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Abstract The behavior of an isolated dilute granular 
gas near the threshold of its clustering instability is in- 
vestigated by means of fluctuating hydrodynamics and 
the direct simulation Monte Carlo method. The theoret- 
ical predictions from the former are shown to be in good 
agreement with the simulation results. The total energy 
of the the system is renormalized by the fluctuations of 
the vorticity field. Moreover, the scaled second moment 
of the energy fluctuations exhibits a power-law divergent 
behavior 

Keywords Homegeous cooling state • clustering 
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| 1 Introduction 

Granular materials often exhibit flows that are similar to 
i those found in molecular fluids. In fact, phenomenolog- 
' ical hydrodynamic-like equations are frequently used to 
describe them [l|. Nevertheless, in general, fluctuations 
of the hydrodynamic fields are much larger in granular 
| systems than in ordinary fluids 0, as a consequence of 
the number of particles being many orders of magnitude 
smaller in the former than in the latter. Thus fluctuations 
are expected to play a more significant role in granular 
flows than in molecular fluid flows. In addition, it must 
, be realized that a new source of noise and fluctuations 
is present in granular systems, associated with the local- 
ized character of the energy dissipation. The properties 
of this intrinsic noise remain largely unknown, although 
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some simple limiting cases have been recently investi- 
gated mil- 

Here, some results for the total energy fluctuations 
in a freely evolving granular gas, below the critical size 
for the onset of the clustering or shear instability [H, 
will be reported. The model considered is a system of TV 
smooth inelastic hard spheres of mass m and diameter a. 
The inelasticity is assumed to be described by a constant, 
velocity independent, coefficient of normal restitution a. 
It is well known that the simplest possible state for this 
system is the so-called homogeneous cooling state (HCS). 
At the level of hydrodynamics, this state is described by 
a constant number density nh, a vanishing flow velocity, 
and a time dependent temperature 2ft (i), that obeys the 
equation 

d t T h (t) = -Ch{t)T h {t), (1) 

where Ch[^h> Tk(t)] is a cooling rate that must be spec- 
ified from a microscopic theory. For the case of hard 
spheres considered here, there is no microscopic energy 
scale associated with the collision model. Therefore, the 
temperature dependence of the cooling rate can be de- 

1/2 

tcrmined by dimensional arguments as Ch(t) oc T h ' (t). 
The HCS is unstable against long wavelength spatial per- 
turbations, leading to the formation of velocity vortices 
and high density clusters of particles [5j . Linear stability 
analysis of the hydrodynamic equations indicates that 
this instability is driven by the transversal shear mode 
[llEl- Using phenomenological Navier-Stokes equations, 
it is found that the critical wavelength L c beyond which 
the system becomes unstable is given by 

*\ 1/2 

- I |-J to- (2) 

; (n^cr 2 ) -1 is proportional to the mean free 
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where Vh(t) = [2T^(t)/m] 1 ^ 2 is a thermal velocity and 
Tf(Th) is the shear viscosity of the granular fluid. For the 
particular case of a dilute gas described by the inelastic 
Boltzmann equation, which is the case to be considered 
in the following, the explicit expressions of r\ and Ch{t) 
are given in [7[. 

Although the initial set up of the clustering instabil- 
ity has been extensively studied, much less attention has 
been payed to the behavior of the system as the insta- 
bility is approached from below, i.e. as the linear system 
size L increases towards L c . The aim of this paper is 
to analyze the behavior of the total energy of of a dilute 
granular gas of hard spheres in that region, by using fluc- 
tuating hydrodynamics and the direct simulation Monte 
Carlo (DSMC) method [![. A similar study for a system 
of hard disks, where the simulations were carried out by 
means of molecular dynamics, has already been reported 
elsewhere |S 03 • 



2 Fluctuating hydrodynamics and average 
energy 

One of the main advantages of the DSMC method [1] 
is that it allows to explote the symmetry of the state 
to be investigated, when dividing the system into cells 
to apply the algorithm. For instance, to study average 
properties of the HCS, the system can be forced to stay 
homogeneous by considering a single cell, the position 
of the particles then becoming irrelevant. This allows 
to significantly increase the statistics of the results, de- 
creasing the uncertainty. Of course, by doing this all the 
spatial hydrodynamic fluctuations and correlations that 
may occur in the system are by definition eliminated. 
This implies, in particular, that the clustering instabil- 
ity can not develop and its influence on the properties of 
the system can not be studied with a single cell. 

Nevertheless, it is possible to consider intermediate 
situations where the clustering instability is present, but 
the system is forced to have some kind of symmetry that 
allows to increase the statistical accuracy. This can be 
achieved, in particular, by dividing the system into a 
given number of parallel layers of the same width. Each 
layer is considered as a cell when applying the DSMC al- 
gorithm. Since the positions of particles inside the same 
cell play no role in determining their collision probability, 
it follows that the dynamics is somehow coarse-grained 
over each layer and variations of the properties inside 
it can not be analyzed. On the other hand, variations 
of the hydrodynamic fields between different layers, i.e. 
along the coordinate perpendicular to them, can show 
up. Therefore, the clustering instability appears when 
the length of the system along that direction, L x , reaches 
the critical value L c . This is the coarse-grained descrip- 
tion for which the theory will be developed in the follow- 
ing. 



The total energy E(t) of the system, either isolated 
or with periodic boundary conditions, we are considering 
can be expressed as 



iV 



mV 2 (t) 



dx 
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^N x (x,t)f{x,t) + ^N x (x,t)u 2 (x,t) 



(5) 

where Vi(t) is the velocity of particle i at time t, the 
tilde is used to indicate that the quantities are treated 
as stochastic variables, and the x axis has been taken 
perpendicular to the layers described above. In the equa- 
tion, N x (x, t) is the number of particles per unit of length 
in the x direction, u(x, t) is the flow velocity field, and 
T(x, t) the temperature. The last two quantities are coarse- 
grained over the cells. The microscopic definitions of 
these fields are 



N x (x,t) 
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N x (x,t)u(x,t)=J2 Vi (t)5 [x - X, (i)] , 



(6) 



(7) 



o N 

-N x (x, t)T(x, t) = ^ ]T [Vi(t) - u(x, t)] 2 5 [x ~ Xi(t)] , 

i=l 

(8) 

that correspond to the idealized limit of infinitely narrow 
layers. Note that, consistently with the discussion above, 
only the x coordinate of the particles at time t, Xi(t), 
appears in the definition of the coarse-grained fields. Ex- 
pansion of Eq. ([5]) around the average values of the fields 
in the HCS, Y a ^ 1 retaining up to second order terms in 
the deviations, SY a (x,t) = Y a ^ — Y a (x,t), yields 

6E{t) = E(t) - E h (t) 

= i j dx {3N x>h 5f(x,t) + 3SN x (x,t)6T(x,t) 

+ mN x ^ h [du(x,t)} 2 }, (9) 

with E h (t) = WT h (t)/2 and N x , h = N/L x . 

In order to eliminate from the analysis the time de- 
pendence of the reference HCS, it is convenient to intro- 
duce dimensionless length £ and time s scales by 



I = 



ds — 



v h (t)dt 



(10) 



V 4) 

Also, dimensionless hydrodynamic fields are defined as 

5N x {x,t) 
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Moreover, the Fourier transform is introduced through 

P*,k{s) = ( dle- m p x {i,s), (14) 



and similarly for the other fields. Then, Eq. ((5]) becomes 
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A x 
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Px,k(s)0-k{s) + — |CCA;(S) | S 



,(15) 



with A x = L x /£q. Then, to calculate e(t), expressions for 
the Fourier components of the fluctuating hydrodynamic 
fields are needed. It will be assumed that, at the meso- 
scopic level used here, they obey Langevin equations ob- 
tained by linearizing the Navier-Stokes equations for a 
granular gas around the HCS. Moreover, if attention is 
restricted to the quasi-elastic limit, i.e. a very close to 
unity, it can be expected that a good approximation be 
obtained by using the same expressions for the properties 
of the noise terms in the Langevin equations as those in 
the Landau-Langevin equations for normal fluids [H|;[l2j]. 

Consider the Fourier transform of the flow field, u>fc(s). 
The vorticity field or transverse flow field in a general 
situation is by definition its component perpendicular to 
the vector fc, i.e. in the present case perpendicular to 
the x axis. It will be denoted by Wfe,j_(s). The coarse- 
grained velocity field uj(£,s) is related with the actual 
velocity field w(l,s) by 



d£ ± w(e,s) = ^-u(e,a) 



(16) 



where I = t/£q, £± denotes the vector component of I 
perpendicular to the x axis, and V* is the volume of the 
system in the dimensionless scale, V* = V/Iq. Using this, 
it is easily obtained that the Landau-Langevin equation 
for LJk,±(s) is 



d_ 

ds~ 



(17) 



The term & ] ± (s) is a Gaussian white noise verifying the 
fluctuation dissipation relation 



A 



N 



(18) 



I being the unit tensor of dimension 2 and the angular 
brackets denoting average over the noise realizations. It 
is worth to remark that the clustering instability mani- 
fests itself in Eq. (|17p. which shows that for \k\ < k c = 

1/2 

(C*/2r7*) ' , the fluctuations of the scaled transverse flow 



field grow in time. On the other hand, for A x < A c 
2w/k c , the long time solution of Eq. (fT7|) is 



J — OO 



(19) 



with Aj_(fc) = ^ - rj*k 2 . Then, by using Eq. (JT5J) it is 
obtained: 

< w fc ,±(s)w fe / I > 



A 2 x y*P 



,(s-s')A x (fc) 



k,-k' 



and, in particular, 
< \u k ±(s)\ 2 >= - 



A 2 xV *k 2 
NX±(k)' 



(20) 



(21) 



Equation ([20]) holds for s > s' > 1. 

A Langevin equation for the energy is obtained by 
linearizing the macroscopic average equation for it, 



--Qh{T h )N x , h I dxST(x,t). 



(22) 



For a normal fluid, the right hand side identically van- 
ishes since the cooling rate is zero. Then, an equation for 
the scaled energy fluctuations is easily found, 



de(s) „ , s 3C* , , 



(23) 



The noise term mentioned above, intrinsic to the local 
energy dissipation in collisions, has been omitted, since it 
is expected to give small contributions as compared with 
those to be kept, in the quasi-elastic limit. Combination 
of Eqs. (HH) and (23|) yields 



tfe(s) 
ds 



= -c 
><E 
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Pk,x(s)0-k(s) 



(24) 



In the threshold of the clustering instability, the fluctu- 
ations of the transversal components of the flow velocity 
are expected to dominate over the density and the longi- 
tudinal velocity fluctuations, so the above equation can 
be reduced to 



rfe(s) 
ds 



(25) 



^From this equation, it follows that the average value of 
the total energy of the system is 



< e 4D<i"fc.x(«)i 2 >-t = -^E v 1 



y z_. Xx (k) 



(26) 



This is the result expected from the expression obtained 
for a system without introducing the coarse-grain aver- 
age over the parallel layers (Toj . 
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For SL = (L c — L x )/L c <C 1 and positive, the main 
contribution to Eq. (|26p is given by those modes hav- 
ing the largest possible wavelength, i.e. those with |fc| = 
kmin — 27r / ' A x . For them, it is 

( 27T ^ 

X±(k) = X±(k min ) = y - rf ( — 



c* 



Lrr: 



-C5L. 



(27) 



It is at this point when the coarse-grain used in the 
simulations deserves some attention. In a normal simu- 
lation, for instance using molecular dynamics of a cubic 
cell, the number of modes with k = k min compatible 
with the periodic boundary conditions should be 6, two 
for each of the axis of the system [l(| . Nevertheless, the 
coarse-grain in the directions perpendicular to the layers, 
actually kills all the modes in those directions. Conse- 
quently only two modes with k — ±k m i n survive. Then, 
from Eqs. {56} and (|27]). 



(28) 



(e) st =A E (5L)-\ A e 



N x , h L c ■ 



where it has been used that near the instability it is 
N ~ N x j t L c . This equation shows that a renormaliza- 
tion by fluctuations of the total energy of the HCS takes 
place as the clustering instability is approached. In fact, 
it is seen that the energy of the system becomes larger 
than the prediction based on Haff's law. When inter- 
preting this theoretical prediction, it must be kept in 
mind the the approximation carried out here only holds 
while the fluctuations are small as compared with their 
average values far from the instability. This requires, in 
particular, that the right hand side of Eq. ([28]) be small. 



3 DSMC method results 

The simulation results to be reported in the following 
have been obtained by using the DSMC method, parti- 
tioning the system in parallel layers as described above. 
The simulations started with a configuration in which 
the particles were homogeneously distributed along the 
length L x of the system and their velocities obeyed a 
Maxwellian distribution. In order to increase the num- 
ber of statistical averages and avoid the technical prob- 
lems inherent to the continuous cooling of the HCS, the 
steady representation of the latter was used. This con- 
sists in a change in the time scale that implies a mod- 
ification of the dynamics, but does not involve any in- 
ternal property of the system. This dynamics leads to a 
steady state whose properties are exactly mapped into 
those of the HCS. The details of the method have been 
extensively analyzed elsewhere [l3| and will be not re- 
peated here. Then, the system is allowed to evolve until 
the steady state is reached. It is important to verify that 
the system has actually reached stationarity, since the 



time required to reach it increases very fast as the in- 
stability is approached. For instance, for a = 0.9, the 
system becomes stationary after about 100 collisions per 
particle for L x /L c = 0.88, while it needs of the order of 
2000 collisions per particle for L x /L c = 0.997. Once the 
system is in the steady state, the different macroscopic 
properties of interest have been computed by taking av- 
erages of the microscopic values. In addition, in all the 
results reported here, the data have also been averaged 
over 1000 independent trajectories. 

The height of the cells used in the simulations has 
been Ax = 0.04£o, significantly smaller than the mean 
free path A = 0.225£o, as required by the DSMC method. 
The number of particles per cell was 100, so N x> h = 
N/L x — 2500£g . The number of cells to be considered 
in each case was determined from the the above figures 
and the value of L x being simulated. It is worth to stress 
that the number of particles does not affect the effective 
dynamics used in the DSMC that remains, by definition, 
that of an iVparticle system in the low density (Boltz- 
mann) limit [8(. 

For a given value of a, the series of simulations were 
started with a system size L x significantly smaller than 
the critical values predicted by Eq. ((2]), namely around 
0.8 times that value. From there, the size of the system 
was increased systematically and its properties studied 
in each case. Special attention was given to verify the ho- 
mogeneity of the hydrodynamic properties of the system. 
Once the formation of velocity vortices was observed, the 
series was stopped since this is the indication that the 
size is larger than the critical length L c . 

In Figs. Q] and [2] the quantity (e)^ 1 is plotted as a 
function of A x — L x /£q for a = 0.9 and a = 0.8, respec- 
tively. In both cases, an increase of the average energy 
relative to the Haff's law prediction is observed as the 
size of the system is decreased approaching the critical 
value. It must be noted, however, that this increase re- 
mains quantitatively small, namely below 10% in all the 
parameter region studied. 

If the theoretical prediction in Eq. |28|) is verified, 
the above plot must lead to a straight line, and that is 
in fact was is obtained when L x increases approaching 
the critical length. This confirms the exponent —1 in Eq. 
([28)) . From the slope and the ordinate at the origin, the 
simulation values of the critical size L c and the critical 
amplitude Ae are obtained. These fitting parameters as 
well as the theoretical predictions for them are given in 
Table [T] The theoretical prediction for the critical length 
has been obtained by means of Eq. ((2]) and the expres- 
sions for the shear viscosity and the cooling rate in the 
first Sonine approximation given in It is observed 
that the agreement between the theory and the simula- 
tion results for the critical size L c is really good. This is 
consistent with the well established result that the first 
Sonine approximation provides good values for the shear 
viscosity and the cooling rate of a dilute granular gas in 
the range of values of a considered here. 
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Fig. 1 Dimensionless quantity {e) st , defined in the main 
text, as a function of the size of the system L x , measured 
in units of £o = (n^cr 2 ) -1 , for a = 0.9. The circles are the 
DSMC results and the solid line a linear fit near the clustering 
instability. 
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Fig. 2 The same as in Fig. [T] but for a system of particles 
with a — 0.8. 



On the other hand, the comparison between the the- 
oretical prediction and the DSMC result for the critical 
amplitude of the energy A e is not so good. In Table Q] 
the measured values for N x L c A e is given. According to 
Eq. (|28p it should be equal to a constant, independent of 
the inelasticity, namely 2. Instead, although the order of 
magnitud is accurately predicted, definitely larger values 
are obtained showing, in addition, a relevant dependence 
on a. Something about the possible origins of this dis- 
crepancy and the improvement of the theory will be said 
at the end of the paper. 



Table 1 Theoretical prediction, and values obtained 

from the simulation data for the average energy, La, and 

from the second moment of the fluctuations, Isf for the crit- 
ical length characterizing the clustering instability of a dilute 
three-dimensional granular g function of the coeffi- 

cient of restitution a. In all cases, L c is measured in units 
of £o = (n/iff 2 )" 1 . The parameters reported in the last two 
columns involve the critical amplitudes A e and A a and should 
be both equal to 2, according to the theoretical predictions. 
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4 Critical energy fluctuations 

Consider next the fluctuations around its average value 
of the total energy of the system, again in the threshold 
of the instability. Since it has been shown above that 
the average value of the energy differs from the Haff law 
prediction, let us define 



S'e(s) = e(s)- < e > st = 

5'uJ(s) = uJ(s) - (p(s)) st 
with 



E(s) - (E(s)) 
E H (s) 



(29) 
(30) 

(31) 



Equation (|26|) shows that < lu > st =< e > s t and, there- 
fore, Eq. ([23)) is equivalent to 

±8>e(s) = -£-[5'e(s)-5>uJ(s)}. (32) 

Taking into account that cumulants of uik j_ of order 

higher than two vanish since £^(s) is assumed to be 
Gaussian, it follows from Eqs. (fl~5|) and (flU)) that 



< 5' uj{s)5' uj(s') >= 



N 2 T 2 



-2(s-s')CSL 



(SLY 



(33) 



Upon deriving this equation, it has been used that only 
modes with k in the direction of the x axis can be excited, 
due to the way in which the simulations are carried out. 
The solution of Eq. (|32|) . once the initial value has been 
forgotten, is 

c* 



5'e(s) = y / ds ie - c * (s - Sl)/2 (5'w(si). 



(34) 



Then, by means of Eq. (|33|) it follows that 



< S'e(s)S'e(s') > = 



D -2(s-s')C 5L 



(SLY 



(35) 
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s > s' S> 1. Terms involving powers of SL larger than —2 
have been consistently neglected in this expression. Thus 
both the total energy fluctuations and the fluctuations of 
the energy associated with the transversal modes decay 
with a characteristic time (2£*<$X) _1 , indicating a diver- 
gent behavior of the relaxation time as SL . Also, the 
second moment of the energy fluctuations is predicted 
to diverge. For s — s' , Eq. ([35]) leads to the stationary 
value 



< {S'ef > st = A 2 JL 



A, = 



V2 



It is 



< (E(a)- < E(s) >) 2 > 
< E(s) > 



N x hL c 



.< {S'ef > 



(36) 



(37) 



To compare with the DSMC method results, it must be 
realized that Eq. (|36|) gives the moment of the energy 
fluctuations that dominate near the instability, but in a 
region where the fluctuations are still small. Then, what 
has been done is to measure u 2 E as a function of L x , start- 
ing from values quite smaller than L c . There, the value 
of a e is not affected by the presence of the instability, 
and Ncr E is independent of L x Q. It will be denoted by 
(Ncr E ) h . Then, the quantity considered has been 



N Xt hL c a E 



-1-1/2 



1 



(38) 



The measured values for this quantity are plotted as a 
function of L. x in Figs. [3] and [H for a — 0.9 and a = 0.8, 
respectively. As predicted by Eq. (|36|) . a linear behav- 
ior is observed as the critical size is approached. Al- 
though the fit is reasonably good for all the values of 
a reported here, it becomes clearly worse as the inelas- 
ticity increases. Moreover, the 'critical' region, identified 
by the values of L x for which the predicted qualitative 
critical behavior is actually observed, is smaller for the 
second moment of the energy than for its average. This 
can be verified by comparing, for instance, Figs. [2] and [4j 
In contrast to what happened with the average energy, 
the growth exhibited by the fluctuations is quite fast. In 
the region plotted in Figs. and |H the second moment 
<je increases more than one order of magnitude. 

^From the slope and ordinate in the origin of the 
fitting straight line, the simulation values of the crit- 
ical length L c and amplitude A a are directly derived. 
The values obtained in this way are also included in Ta- 
ble [TJ Again, the results for the critical length L c are 
in excellent agreement with the theoretical predictions. 
Moreover, the fact that the same values follow from the 
measurements of both the average energy and the second 
moments provides a test of the internal consistency of the 
theory. With regard to the critical amplitude A a , signif- 
icant deviations from the the theoretical predictions are 
also found in this case, although they are smaller than 
for A e . In fact, for the most elastic system considered, 
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Fig. 3 Plot of the DSMC re sult s (circles) for the dimension- 
less quantity denned in Eq. (|38[1 as a function of the size of 
the system L x , for a — 0.9. The solid line is a linear fit near 
the clustering instability. 



a = 0.95, there is a good agreement. The most relevant 
disagreement between theory and simulations is the def- 
inite dependence on inelasticity of the parameter combi- 
nation given in Table [T] shown by the latter, while the 
former predicts a constant value (namely, 2). 

It can be wondered which is the aspect of the theoret- 
ical approach developed here that should be modified in 
order to improve the accuracy of the predicted expres- 
sions for the critical amplitudes. Of course, a first im- 
portant limitation of the theory is its restriction to the 
almost elastic limit, implied by the use of the Landau 
fluctuation-dissipation relations. But it must be realized 
that up to now those expressions have not been gener- 
alized for non-conservative interactions. A more modest 
step in this direction should be to include in Eq. (|22|) the 
intrinsic noise term associated with the cooling rate, as 
already indicated. 

In refs. @ and [l(| , a scaling property of the probabil- 
ity distribution of the energy fluctuations in the thresh- 
old of the clustering instability of a two-dimensional gran- 
ular fluid was identified. Moreover, the scaling function 
was very well fitted by the same expression as found in 
several equilibrium and non-equilibrium molecular sys- 
tems [HI EH ■ The possible existence of a similar scaling 
for the three-dimensional granular gas considered here 
has also been investigated, finding similar results. 
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